Version 1.0, July 2025, SA
This script calculates a weighted average of data clustered together 100% consistently across leave1out iterations, and runs the statistical comparisons for each supercluster.
It plots the supercluster connectivity profiles in Figure 4.
It also runs the statistical analyses and plots for inside vs outside infantile amnesia window comparisons (Figure 5)
Input: threshMat.csv and Hippo_BinxAxis_F_betas.txt
library(dplyr)
library(stringr)
library(emmeans)
library(ggplot2)
library(nlme)
library(car)
library(Rmisc)
### specify directory, CHANGE TO YOUR PATH
dir<-"/Users/audrainsp/Library/CloudStorage/OneDrive-NationalInstitutesofHealth/BabyHippos/R_BabyHippos"
#### read in binarized, thresholded matrix.
### 1 is where clustering was consistent across all subjects, 0 was inconsistent
data_threshMat <- read.csv(paste(dir, "/threshMat.csv", sep=""), header = FALSE)
### get rid of diagonal and lower triangle of the binarized matrix
diag(data_threshMat)<-0
data_threshMat[lower.tri(data_threshMat)]<-0
### read in our original data and format it for our purposes
data <- read.table(paste(dir, "/Hippo_BinxAxis_F_betas.txt", sep=""), header = TRUE)
### keep the top 44 clusters for cluster threshold of 10
data<-subset(data, clust < 45)
## collapse over hem for each subj
data_ag_hem<-aggregate(data[, "beta"], by=(list(data$subj, data$ax, data$bin, data$clust)), mean)
colnames(data_ag_hem)<-c("subj","ax","bin","clust","beta")
#### specify which clusters/rois go into each supercluster
supercluster_1<-c(1,6,7,28,34,36,41,43)
supercluster_2<-c(2,3,9,11,13,14,16,19,20,22,24,30,31,32,37,38,40,44)
supercluster_3<-c(4,15,17,23,25,26,29,33,35,39)
supercluster_4<-c(5,10)
supercluster_5<-c(8,12,27,42)
supercluster_6<-c(18,21)
### append supercluster info to our dataframe
data_ag_hem<-mutate(data_ag_hem, supercluster = case_when(clust %in% supercluster_1 ~ 1,
clust %in% supercluster_2 ~ 2,
clust %in% supercluster_3 ~ 3,
clust %in% supercluster_4 ~ 4,
clust %in% supercluster_5 ~ 5,
clust %in% supercluster_6 ~ 6))
### loop through to create a column of ages and append to our dataframe
ages<-as.data.frame(c("00mo", "01mo", "02mo", "03mo", "04mo", "05mo", "06mo", "07mo", "08mo", "09mo", "10mo", "11mo", "12mo", "13mo", "14mo", "15mo", "16mo", "17mo", "18mo", "19mo", "20mo", "21mo", "22mo", "23mo", "24mo", "25mo", "26mo"))
colnames(ages)<-c("age")
for (r in 1:nrow(data_ag_hem)){
for (a in 1:nrow(ages))
if ((str_detect(data_ag_hem[r,"subj"], pattern = ages[a,])) ==TRUE) {
data_ag_hem[r,"age"]<-a-1
} else if ((str_detect(data_ag_hem[r,"subj"], pattern = "00mo"))) {
data_ag_hem[r,"age"]<-0
}
}
### initialize some empty lists to append to later, for each supercluster
data_sc_1<-data.frame(matrix(nrow=0, ncol=6))
colnames(data_sc_1)<-c("subj","bin","age","ax","superclust","beta")
data_sc_2<-data.frame(matrix(nrow=0, ncol=6))
colnames(data_sc_2)<-c("subj","bin","age","ax","superclust","beta")
data_sc_3<-data.frame(matrix(nrow=0, ncol=6))
colnames(data_sc_3)<-c("subj","bin","age","ax","superclust","beta")
data_sc_4<-data.frame(matrix(nrow=0, ncol=6))
colnames(data_sc_4)<-c("subj","bin","age","ax","superclust","beta")
data_sc_5<-data.frame(matrix(nrow=0, ncol=6))
colnames(data_sc_5)<-c("subj","bin","age","ax","superclust","beta")
data_sc_6<-data.frame(matrix(nrow=0, ncol=6))
colnames(data_sc_6)<-c("subj","bin","age","ax","superclust","beta")
### loop through matrix and when there is a 1, average the two roi hippo-cortical conn values together
for (r in 1:nrow(data_threshMat)) {
for (c in 1:ncol(data_threshMat)) {
if (data_threshMat[r,c]==1){
#print("yes")
temp_subset<-subset(data_ag_hem, (clust==r | clust==c))
### average across clust
temp_agg<-aggregate(temp_subset[, "beta"], by=(list(temp_subset$subj,temp_subset$bin, temp_subset$age, temp_subset$ax, temp_subset$supercluster)), mean)
colnames(temp_agg)<-c("subj","bin","age","ax","superclust","beta")
if (temp_agg[1,"superclust"]=="1") {
data_sc_1<-rbind(data_sc_1, temp_agg)
}
if (temp_agg[1,"superclust"]=="2") {
data_sc_2<-rbind(data_sc_2, temp_agg)
}
if (temp_agg[1,"superclust"]=="3") {
data_sc_3<-rbind(data_sc_3, temp_agg)
}
if (temp_agg[1,"superclust"]=="4") {
data_sc_4<-rbind(data_sc_4, temp_agg)
}
if (temp_agg[1,"superclust"]=="5") {
data_sc_5<-rbind(data_sc_5, temp_agg)
}
if (temp_agg[1,"superclust"]=="6") {
data_sc_6<-rbind(data_sc_6, temp_agg)
}
}
}
}
### data for each supercluster weighted by consistency across leave-one-out analysis
data_sc_1<-aggregate(data_sc_1[, "beta"], by=(list(data_sc_1$subj, data_sc_1$bin, data_sc_1$age, data_sc_1$ax)), mean)
colnames(data_sc_1)<-c("subj","bin","age","ax","beta")
data_sc_1$supercluster<-1
data_sc_2<-aggregate(data_sc_2[, "beta"], by=(list(data_sc_2$subj, data_sc_2$bin, data_sc_2$age, data_sc_2$ax)), mean)
colnames(data_sc_2)<-c("subj","bin","age","ax","beta")
data_sc_2$supercluster<-2
data_sc_3<-aggregate(data_sc_3[, "beta"], by=(list(data_sc_3$subj, data_sc_3$bin, data_sc_3$age, data_sc_3$ax)), mean)
colnames(data_sc_3)<-c("subj","bin","age","ax","beta")
data_sc_3$supercluster<-3
## nothing in this one, parahippocampal. Because the ROIs were not clustered together 100% of the time.
# data_sc_4<-aggregate(data_sc_4[, "beta"], by=(list(data_sc_4$subj, data_sc_4$bin)), mean)
# colnames(data_sc_4)<-c("subj","bin","beta")
# data_sc_4$supercluster<-4
data_sc_5<-aggregate(data_sc_5[, "beta"], by=(list(data_sc_5$subj, data_sc_5$bin, data_sc_5$age, data_sc_5$ax)), mean)
colnames(data_sc_5)<-c("subj","bin","age","ax","beta")
data_sc_5$supercluster<-5
data_sc_6<-aggregate(data_sc_6[, "beta"], by=(list(data_sc_6$subj, data_sc_6$bin, data_sc_6$age, data_sc_6$ax)), mean)
colnames(data_sc_6)<-c("subj","bin","age","ax","beta")
data_sc_6$supercluster<-6
#data_LOO<-rbind(data_sc_1,data_sc_2,data_sc_3,data_sc_5,data_sc_6)
data_SC_all<-rbind(data_sc_1, data_sc_2, data_sc_3, data_sc_5, data_sc_6)
### turn variables into factors
data_SC_all$bin<-as.factor(data_SC_all$bin)
data_SC_all$supercluster<-as.factor(data_SC_all$supercluster)
## model
## this model takes forever to run, so can also just read a saved version in using readRDS
## model is weighted to allow heteroskedasticity
#model<-lme(beta ~ ax*bin*supercluster, random= ~1|subj, weights=varIdent(form=~1|bin*ax*supercluster), data=data_SC_all, control = lmeControl(msMaxIter=1000, msMaxEval=1000), na.action=na.omit)
model<-readRDS('/Users/audrainsp/Library/CloudStorage/OneDrive-NationalInstitutesofHealth/BabyHippos/R_BabyHippos/GitHub/mod1_weighted_ax_bin_supercluster.rda') ### CHANGE TO YOUR PATH. OR RUN ABOVE MODEL.
## model diagnostics
hist(residuals(model))
plot(model)
## stats
anova(model)
#### post hocs
## these are the pairwise comparisons used in the manuscript
## stats for 2 way interactions for each supercluster can be found next to each subplot below
emmeans(model, list(pairwise ~ ax|bin|supercluster), adjust = "none")
$`emmeans of ax | bin, supercluster`
bin = 1, supercluster = 1:
ax emmean SE df lower.CL upper.CL
ant 0.0667 0.00486 211 0.0572 0.0763
post 0.0379 0.00500 211 0.0280 0.0477
bin = 2, supercluster = 1:
ax emmean SE df lower.CL upper.CL
ant 0.0642 0.00393 208 0.0564 0.0719
post 0.0543 0.00406 208 0.0463 0.0623
bin = 3, supercluster = 1:
ax emmean SE df lower.CL upper.CL
ant 0.0660 0.00416 208 0.0578 0.0742
post 0.0623 0.00407 208 0.0543 0.0703
bin = 4, supercluster = 1:
ax emmean SE df lower.CL upper.CL
ant 0.0458 0.00457 208 0.0368 0.0549
post 0.0599 0.00450 208 0.0510 0.0687
bin = 1, supercluster = 2:
ax emmean SE df lower.CL upper.CL
ant 0.0432 0.00443 211 0.0345 0.0519
post 0.0516 0.00482 211 0.0421 0.0611
bin = 2, supercluster = 2:
ax emmean SE df lower.CL upper.CL
ant 0.0332 0.00499 208 0.0233 0.0430
post 0.0736 0.00494 208 0.0639 0.0834
bin = 3, supercluster = 2:
ax emmean SE df lower.CL upper.CL
ant 0.0298 0.00554 208 0.0189 0.0408
post 0.0602 0.00491 208 0.0506 0.0699
bin = 4, supercluster = 2:
ax emmean SE df lower.CL upper.CL
ant 0.0335 0.00494 208 0.0238 0.0433
post 0.0590 0.00481 208 0.0495 0.0685
bin = 1, supercluster = 3:
ax emmean SE df lower.CL upper.CL
ant 0.0697 0.00542 211 0.0590 0.0804
post 0.0778 0.00528 211 0.0674 0.0882
bin = 2, supercluster = 3:
ax emmean SE df lower.CL upper.CL
ant 0.1113 0.00517 208 0.1011 0.1214
post 0.0818 0.00488 208 0.0722 0.0914
bin = 3, supercluster = 3:
ax emmean SE df lower.CL upper.CL
ant 0.1098 0.00548 208 0.0990 0.1206
post 0.0851 0.00550 208 0.0743 0.0960
bin = 4, supercluster = 3:
ax emmean SE df lower.CL upper.CL
ant 0.0906 0.00579 208 0.0792 0.1020
post 0.0808 0.00545 208 0.0700 0.0915
bin = 1, supercluster = 5:
ax emmean SE df lower.CL upper.CL
ant 0.0637 0.00615 211 0.0516 0.0758
post 0.1031 0.00707 211 0.0892 0.1171
bin = 2, supercluster = 5:
ax emmean SE df lower.CL upper.CL
ant 0.1132 0.00588 208 0.1016 0.1248
post 0.1196 0.00572 208 0.1083 0.1309
bin = 3, supercluster = 5:
ax emmean SE df lower.CL upper.CL
ant 0.1054 0.00492 208 0.0957 0.1151
post 0.1066 0.00529 208 0.0962 0.1171
bin = 4, supercluster = 5:
ax emmean SE df lower.CL upper.CL
ant 0.0953 0.00575 208 0.0839 0.1066
post 0.1047 0.00605 208 0.0928 0.1166
bin = 1, supercluster = 6:
ax emmean SE df lower.CL upper.CL
ant 0.2451 0.00915 211 0.2271 0.2631
post 0.1383 0.00657 211 0.1253 0.1512
bin = 2, supercluster = 6:
ax emmean SE df lower.CL upper.CL
ant 0.2449 0.00735 208 0.2304 0.2594
post 0.1731 0.00750 208 0.1583 0.1878
bin = 3, supercluster = 6:
ax emmean SE df lower.CL upper.CL
ant 0.2294 0.00737 208 0.2149 0.2439
post 0.1647 0.00687 208 0.1511 0.1782
bin = 4, supercluster = 6:
ax emmean SE df lower.CL upper.CL
ant 0.1936 0.00909 208 0.1757 0.2115
post 0.1386 0.00725 208 0.1243 0.1529
Degrees-of-freedom method: containment
Confidence level used: 0.95
$`pairwise differences of ax | bin, supercluster`
bin = 1, supercluster = 1:
3 estimate SE df t.ratio p.value
ant - post 0.02888 0.00437 1872 6.603 <.0001
bin = 2, supercluster = 1:
3 estimate SE df t.ratio p.value
ant - post 0.00987 0.00263 1872 3.753 0.0002
bin = 3, supercluster = 1:
3 estimate SE df t.ratio p.value
ant - post 0.00372 0.00291 1872 1.280 0.2006
bin = 4, supercluster = 1:
3 estimate SE df t.ratio p.value
ant - post -0.01400 0.00288 1872 -4.858 <.0001
bin = 1, supercluster = 2:
3 estimate SE df t.ratio p.value
ant - post -0.00846 0.00365 1872 -2.321 0.0204
bin = 2, supercluster = 2:
3 estimate SE df t.ratio p.value
ant - post -0.04049 0.00493 1872 -8.214 <.0001
bin = 3, supercluster = 2:
3 estimate SE df t.ratio p.value
ant - post -0.03039 0.00541 1872 -5.617 <.0001
bin = 4, supercluster = 2:
3 estimate SE df t.ratio p.value
ant - post -0.02549 0.00383 1872 -6.647 <.0001
bin = 1, supercluster = 3:
3 estimate SE df t.ratio p.value
ant - post -0.00812 0.00527 1872 -1.541 0.1235
bin = 2, supercluster = 3:
3 estimate SE df t.ratio p.value
ant - post 0.02946 0.00505 1872 5.828 <.0001
bin = 3, supercluster = 3:
3 estimate SE df t.ratio p.value
ant - post 0.02469 0.00590 1872 4.183 <.0001
bin = 4, supercluster = 3:
3 estimate SE df t.ratio p.value
ant - post 0.00979 0.00552 1872 1.775 0.0761
bin = 1, supercluster = 5:
3 estimate SE df t.ratio p.value
ant - post -0.03944 0.00763 1872 -5.168 <.0001
bin = 2, supercluster = 5:
3 estimate SE df t.ratio p.value
ant - post -0.00639 0.00650 1872 -0.983 0.3257
bin = 3, supercluster = 5:
3 estimate SE df t.ratio p.value
ant - post -0.00128 0.00517 1872 -0.248 0.8039
bin = 4, supercluster = 5:
3 estimate SE df t.ratio p.value
ant - post -0.00942 0.00608 1872 -1.550 0.1213
bin = 1, supercluster = 6:
3 estimate SE df t.ratio p.value
ant - post 0.10681 0.00987 1872 10.819 <.0001
bin = 2, supercluster = 6:
3 estimate SE df t.ratio p.value
ant - post 0.07187 0.00923 1872 7.785 <.0001
bin = 3, supercluster = 6:
3 estimate SE df t.ratio p.value
ant - post 0.06473 0.00872 1872 7.422 <.0001
bin = 4, supercluster = 6:
3 estimate SE df t.ratio p.value
ant - post 0.05500 0.01010 1872 5.435 <.0001
Degrees-of-freedom method: containment
emmeans(model, list(pairwise ~ bin|ax|supercluster), adjust = "none")
$`emmeans of bin | ax, supercluster`
ax = ant, supercluster = 1:
bin emmean SE df lower.CL upper.CL
1 0.0667 0.00486 211 0.0572 0.0763
2 0.0642 0.00393 208 0.0564 0.0719
3 0.0660 0.00416 208 0.0578 0.0742
4 0.0458 0.00457 208 0.0368 0.0549
ax = post, supercluster = 1:
bin emmean SE df lower.CL upper.CL
1 0.0379 0.00500 211 0.0280 0.0477
2 0.0543 0.00406 208 0.0463 0.0623
3 0.0623 0.00407 208 0.0543 0.0703
4 0.0599 0.00450 208 0.0510 0.0687
ax = ant, supercluster = 2:
bin emmean SE df lower.CL upper.CL
1 0.0432 0.00443 211 0.0345 0.0519
2 0.0332 0.00499 208 0.0233 0.0430
3 0.0298 0.00554 208 0.0189 0.0408
4 0.0335 0.00494 208 0.0238 0.0433
ax = post, supercluster = 2:
bin emmean SE df lower.CL upper.CL
1 0.0516 0.00482 211 0.0421 0.0611
2 0.0736 0.00494 208 0.0639 0.0834
3 0.0602 0.00491 208 0.0506 0.0699
4 0.0590 0.00481 208 0.0495 0.0685
ax = ant, supercluster = 3:
bin emmean SE df lower.CL upper.CL
1 0.0697 0.00542 211 0.0590 0.0804
2 0.1113 0.00517 208 0.1011 0.1214
3 0.1098 0.00548 208 0.0990 0.1206
4 0.0906 0.00579 208 0.0792 0.1020
ax = post, supercluster = 3:
bin emmean SE df lower.CL upper.CL
1 0.0778 0.00528 211 0.0674 0.0882
2 0.0818 0.00488 208 0.0722 0.0914
3 0.0851 0.00550 208 0.0743 0.0960
4 0.0808 0.00545 208 0.0700 0.0915
ax = ant, supercluster = 5:
bin emmean SE df lower.CL upper.CL
1 0.0637 0.00615 211 0.0516 0.0758
2 0.1132 0.00588 208 0.1016 0.1248
3 0.1054 0.00492 208 0.0957 0.1151
4 0.0953 0.00575 208 0.0839 0.1066
ax = post, supercluster = 5:
bin emmean SE df lower.CL upper.CL
1 0.1031 0.00707 211 0.0892 0.1171
2 0.1196 0.00572 208 0.1083 0.1309
3 0.1066 0.00529 208 0.0962 0.1171
4 0.1047 0.00605 208 0.0928 0.1166
ax = ant, supercluster = 6:
bin emmean SE df lower.CL upper.CL
1 0.2451 0.00915 211 0.2271 0.2631
2 0.2449 0.00735 208 0.2304 0.2594
3 0.2294 0.00737 208 0.2149 0.2439
4 0.1936 0.00909 208 0.1757 0.2115
ax = post, supercluster = 6:
bin emmean SE df lower.CL upper.CL
1 0.1383 0.00657 211 0.1253 0.1512
2 0.1731 0.00750 208 0.1583 0.1878
3 0.1647 0.00687 208 0.1511 0.1782
4 0.1386 0.00725 208 0.1243 0.1529
Degrees-of-freedom method: containment
Confidence level used: 0.95
$`pairwise differences of bin | ax, supercluster`
ax = ant, supercluster = 1:
3 estimate SE df t.ratio p.value
bin1 - bin2 0.002574 0.00625 208 0.412 0.6809
bin1 - bin3 0.000707 0.00640 208 0.111 0.9121
bin1 - bin4 0.020894 0.00667 208 3.133 0.0020
bin2 - bin3 -0.001867 0.00573 208 -0.326 0.7447
bin2 - bin4 0.018320 0.00603 208 3.039 0.0027
bin3 - bin4 0.020187 0.00618 208 3.267 0.0013
ax = post, supercluster = 1:
3 estimate SE df t.ratio p.value
bin1 - bin2 -0.016433 0.00644 208 -2.551 0.0115
bin1 - bin3 -0.024447 0.00645 208 -3.789 0.0002
bin1 - bin4 -0.021989 0.00673 208 -3.268 0.0013
bin2 - bin3 -0.008014 0.00575 208 -1.394 0.1648
bin2 - bin4 -0.005557 0.00606 208 -0.917 0.3600
bin3 - bin4 0.002458 0.00607 208 0.405 0.6859
ax = ant, supercluster = 2:
3 estimate SE df t.ratio p.value
bin1 - bin2 0.010028 0.00667 208 1.504 0.1341
bin1 - bin3 0.013335 0.00709 208 1.881 0.0613
bin1 - bin4 0.009667 0.00663 208 1.458 0.1464
bin2 - bin3 0.003306 0.00745 208 0.444 0.6577
bin2 - bin4 -0.000361 0.00702 208 -0.052 0.9590
bin3 - bin4 -0.003668 0.00742 208 -0.494 0.6216
ax = post, supercluster = 2:
3 estimate SE df t.ratio p.value
bin1 - bin2 -0.021998 0.00690 208 -3.187 0.0017
bin1 - bin3 -0.008597 0.00688 208 -1.250 0.2126
bin1 - bin4 -0.007362 0.00681 208 -1.082 0.2806
bin2 - bin3 0.013401 0.00696 208 1.924 0.0557
bin2 - bin4 0.014636 0.00690 208 2.122 0.0350
bin3 - bin4 0.001235 0.00687 208 0.180 0.8575
ax = ant, supercluster = 3:
3 estimate SE df t.ratio p.value
bin1 - bin2 -0.041533 0.00749 208 -5.546 <.0001
bin1 - bin3 -0.040109 0.00771 208 -5.205 <.0001
bin1 - bin4 -0.020849 0.00793 208 -2.630 0.0092
bin2 - bin3 0.001425 0.00753 208 0.189 0.8501
bin2 - bin4 0.020685 0.00776 208 2.666 0.0083
bin3 - bin4 0.019260 0.00797 208 2.417 0.0165
ax = post, supercluster = 3:
3 estimate SE df t.ratio p.value
bin1 - bin2 -0.003959 0.00719 208 -0.551 0.5826
bin1 - bin3 -0.007300 0.00763 208 -0.957 0.3395
bin1 - bin4 -0.002946 0.00759 208 -0.388 0.6983
bin2 - bin3 -0.003341 0.00736 208 -0.454 0.6502
bin2 - bin4 0.001013 0.00732 208 0.138 0.8901
bin3 - bin4 0.004354 0.00775 208 0.562 0.5747
ax = ant, supercluster = 5:
3 estimate SE df t.ratio p.value
bin1 - bin2 -0.049474 0.00851 208 -5.817 <.0001
bin1 - bin3 -0.041650 0.00788 208 -5.289 <.0001
bin1 - bin4 -0.031583 0.00842 208 -3.752 0.0002
bin2 - bin3 0.007824 0.00767 208 1.021 0.3086
bin2 - bin4 0.017890 0.00822 208 2.176 0.0307
bin3 - bin4 0.010067 0.00757 208 1.330 0.1850
ax = post, supercluster = 5:
3 estimate SE df t.ratio p.value
bin1 - bin2 -0.016425 0.00909 208 -1.806 0.0724
bin1 - bin3 -0.003492 0.00883 208 -0.396 0.6928
bin1 - bin4 -0.001561 0.00931 208 -0.168 0.8670
bin2 - bin3 0.012932 0.00779 208 1.659 0.0985
bin2 - bin4 0.014864 0.00833 208 1.784 0.0759
bin3 - bin4 0.001931 0.00804 208 0.240 0.8104
ax = ant, supercluster = 6:
3 estimate SE df t.ratio p.value
bin1 - bin2 0.000161 0.01170 208 0.014 0.9891
bin1 - bin3 0.015677 0.01170 208 1.334 0.1836
bin1 - bin4 0.051511 0.01290 208 3.994 0.0001
bin2 - bin3 0.015516 0.01040 208 1.491 0.1375
bin2 - bin4 0.051350 0.01170 208 4.394 <.0001
bin3 - bin4 0.035835 0.01170 208 3.063 0.0025
ax = post, supercluster = 6:
3 estimate SE df t.ratio p.value
bin1 - bin2 -0.034779 0.00997 208 -3.488 0.0006
bin1 - bin3 -0.026401 0.00951 208 -2.776 0.0060
bin1 - bin4 -0.000294 0.00979 208 -0.030 0.9761
bin2 - bin3 0.008378 0.01020 208 0.824 0.4110
bin2 - bin4 0.034486 0.01040 208 3.305 0.0011
bin3 - bin4 0.026108 0.00999 208 2.613 0.0096
Degrees-of-freedom method: containment
### save the output, if desired
# capture.output(emmeans(model, list(pairwise ~ bin|ax|supercluster), adjust = "fdr"), file="/Users/samanthaaudrain/Documents/Science_Projects/BabyHippos/R_BabyHippos/mod2_weighted_ax_bin_supercluster_emmeans1_fdr.doc")
#
# capture.output(emmeans(model, list(pairwise ~ ax|bin|supercluster), adjust = "fdr"), file="/Users/samanthaaudrain/Documents/Science_Projects/BabyHippos/R_BabyHippos/mod2_weighted_ax_bin_supercluster_emmeans2_fdr.doc")
## prepare the data
data_SC<-data_sc_1 ### this is weighted averaged data across component clusters
data_SC$bin<-as.factor(data_SC$bin)
## run the model
model<-lme(beta ~ ax*bin, random= ~1|subj, weights=varIdent(form=~1|bin), data=data_SC, na.action=na.omit)
## model diagnostics
# hist(residuals(model))
# plot(model)
# stats
anova(model)
## summarize for plotting
stats<-summarySE(data=data_SC, measurevar = "beta", groupvars = c("bin","ax"), na.rm=TRUE)
stats$upper_se<-stats$beta + stats$se
stats$lower_se<-stats$beta - stats$se
## plot
#tiff("~/Downloads/cing_operc.tiff", units="in", width=4.5, height=4, res=1000)
ggplot(stats, aes(x=bin, y=beta, group=ax)) +
geom_line(size=1, aes(linetype="solid", color=ax)) +
geom_errorbar(width = .0, aes(ymin=lower_se, ymax=upper_se)) +
theme_classic() +
theme(axis.title = element_text(size=24), axis.ticks.x=element_blank(), text = element_text(size = 18), axis.text.x = element_text(size=20), axis.text.y = element_text(size=22), plot.title = element_text(face="bold", hjust=0.5, size=22), legend.position="none")+
scale_x_discrete(limits=c("1", "2", "3", "4"), labels=c("0-6", "7-12", "13-18","19-25"))+
scale_y_continuous(expand = c(0, 0), limits = c(0,0.3), breaks=seq(0,0.3,.05))+
scale_color_manual(values=c('#152370','#6577a1')) +
geom_hline(yintercept=0)+
# ggtitle("cingulo-opercular
# supercluster")+
labs(x="age", y="connectivity", color="hippocampal axis")
#dev.off()
#ggsave("cingulo-operc.png", path="~/Downloads", width = 4.5, height = 4)
## prepare the data
data_SC<-data_sc_2 ### this is weighted averaged data across component clusters
data_SC$bin<-as.factor(data_SC$bin)
## run the model
model<-lme(beta ~ ax*bin, random= ~1|subj, data=data_SC, na.action=na.omit)
## model diagnostics
# hist(residuals(model))
# plot(model)
## stats
anova(model)
## summarize for plotting
stats<-summarySE(data=data_SC, measurevar = "beta", groupvars = c("bin","ax"), na.rm=TRUE)
stats$upper_se<-stats$beta + stats$se
stats$lower_se<-stats$beta - stats$se
## plot
#tiff("~/Downloads/dFPN.tiff", units="in", width=4.5, height=4, res=1000)
ggplot(stats, aes(x=bin, y=beta, group=ax)) +
geom_line(size=1, aes(linetype="solid", color=ax)) +
geom_errorbar(width = .0, aes(ymin=lower_se, ymax=upper_se)) +
theme_classic() +
theme(axis.title = element_text(size=24), axis.ticks.x=element_blank(), text = element_text(size = 18), axis.text.x = element_text(size=20), axis.text.y = element_text(size=22), plot.title = element_text(face="bold", hjust=0.5, size=22), legend.position="none")+
scale_x_discrete(limits=c("1", "2", "3", "4"), labels=c("0-6", "7-12", "13-18","19-25"))+
scale_y_continuous(expand = c(0, 0), limits = c(0,0.3), breaks=seq(0,0.3,.05))+
scale_color_manual(values=c('#4bf52a','#bdf7b2')) +
geom_hline(yintercept=0)+
# ggtitle("dorsal frontal parietal
# supercluster")+
labs(x="age", y="connectivity", color="hippocampal axis")
#dev.off()
#ggsave("dorsalFPN.png", path="~/Downloads", width = 4.5, height = 4)
## prepare the data
data_SC<-data_sc_3 ### this is weighted averaged data across component clusters
data_SC$bin<-as.factor(data_SC$bin)
## run the model
model<-lme(beta ~ ax*bin, random= ~1|subj, data=data_SC, na.action=na.omit)
## model diagnostics
# hist(residuals(model))
# plot(model)
## stats
anova(model)
## summarize for plotting
stats<-summarySE(data=data_SC, measurevar = "beta", groupvars = c("bin","ax"), na.rm=TRUE)
stats$upper_se<-stats$beta + stats$se
stats$lower_se<-stats$beta - stats$se
## plot
#tiff("~/Downloads/mPFCSTS.tiff", units="in", width=4.5, height=4, res=1000)
ggplot(stats, aes(x=bin, y=beta, group=ax)) +
geom_line(size=1, aes(linetype="solid", color=ax)) +
geom_errorbar(width = .0, aes(ymin=lower_se, ymax=upper_se)) +
theme_classic() +
theme(axis.title = element_text(size=24), axis.ticks.x=element_blank(), text = element_text(size = 18), axis.text.x = element_text(size=20), axis.text.y = element_text(size=22), plot.title = element_text(face="bold", hjust=0.5, size=22), legend.position="none")+
scale_x_discrete(limits=c("1", "2", "3", "4"), labels=c("0-6", "7-12", "13-18","19-25"))+
scale_y_continuous(expand = c(0, 0), limits = c(0,0.3), breaks=seq(0,0.3,.05))+
scale_color_manual(values=c('#f2f207','#fafab4')) +
geom_hline(yintercept=0)+
# ggtitle("vmPFC-STS
# supercluster")+
labs(x="age", y="connectivity", color="hippocampal axis")
#dev.off()
#ggsave("vmpfc_sts.png", path="~/Downloads", width = 4.5, height = 4)
## prepare the data
data_SC<-data_sc_5 ### this is weighted averaged data across component clusters
data_SC$bin<-as.factor(data_SC$bin)
## run the model
model<-lme(beta ~ ax*bin, random= ~1|subj, data=data_SC, na.action=na.omit)
## model diagnostics
# hist(residuals(model))
# plot(model)
## stats
anova(model)
## summarize for plotting
stats<-summarySE(data=data_SC, measurevar = "beta", groupvars = c("bin","ax"), na.rm=TRUE)
stats$upper_se<-stats$beta + stats$se
stats$lower_se<-stats$beta - stats$se
## plot
#tiff("~/Downloads/parietal.tiff", units="in", width=4.5, height=4, res=1000)
ggplot(stats, aes(x=bin, y=beta, group=ax)) +
geom_line(size=1, aes(linetype="solid", color=ax)) +
geom_errorbar(width = .0, aes(ymin=lower_se, ymax=upper_se)) +
theme_classic() +
theme(axis.title = element_text(size=24), axis.ticks.x=element_blank(), text = element_text(size = 18), axis.text.x = element_text(size=20), axis.text.y = element_text(size=22), plot.title = element_text(face="bold", hjust=0.5, size=22), legend.position="none")+
scale_x_discrete(limits=c("1", "2", "3", "4"), labels=c("0-6", "7-12", "13-18","19-25"))+
scale_y_continuous(expand = c(0, 0), limits = c(0,0.3), breaks=seq(0,0.3,.05))+
scale_color_manual(values=c('#fa750f','#faa666')) +
geom_hline(yintercept=0)+
# ggtitle("medial parietal
# supercluster")+
labs(x="age", y="connectivity", color="hippocampal axis")
#dev.off()
#ggsave("medialparietal.png", path="~/Downloads", width = 4.5, height = 4)
## prepare the data
data_SC<-data_sc_6 ### this is weighted averaged data across component clusters
data_SC$bin<-as.factor(data_SC$bin)
## run the model
model<-lme(beta ~ ax*bin, random= ~1|subj, data=data_SC, na.action=na.omit)
## model diagnostics
# hist(residuals(model))
# plot(model)
## stats
anova(model)
## summarize for plotting
stats<-summarySE(data=data_SC, measurevar = "beta", groupvars = c("bin","ax"), na.rm=TRUE)
stats$upper_se<-stats$beta + stats$se
stats$lower_se<-stats$beta - stats$se
## plot
#tiff("~/Downloads/EC.tiff", units="in", width=4.5, height=4, res=1000)
ggplot(stats, aes(x=bin, y=beta, group=ax)) +
geom_line(size=1, aes(linetype="solid", color=ax)) +
geom_errorbar(width = .0, aes(ymin=lower_se, ymax=upper_se)) +
theme_classic() +
theme(axis.title = element_text(size=24), axis.ticks.x=element_blank(), text = element_text(size = 18), axis.text.x = element_text(size=20), axis.text.y = element_text(size=22), plot.title = element_text(face="bold", hjust=0.5, size=22), legend.position="none")+
scale_x_discrete(limits=c("1", "2", "3", "4"), labels=c("0-6", "7-12", "13-18","19-25"))+
scale_y_continuous(expand = c(0, 0), limits = c(0,0.3), breaks=seq(0,0.3,.05))+
scale_color_manual(values=c('#28fcfc','#bfffff')) +
geom_hline(yintercept=0)+
# ggtitle("entorhinal
# supercluster")+
labs(x="age", y="connectivity", color="hippocampal axis")
#dev.off()
#ggsave("entorhinal.png", path="~/Downloads", width = 4.5, height = 4)
Plotting anterior and posterior data separately and overlaid for best visualization
SC_data<-subset(data_sc_1, supercluster==1 & ax=="ant")
#tiff("~/Downloads/CP_1.tiff", units="in", width=4.5, height=4, res=1000)
ggplot(SC_data, aes(y=beta, x=age)) +
geom_vline(xintercept = 6, color="grey", linetype="solid") +
geom_vline(xintercept = 12, color="grey", linetype="solid") +
geom_vline(xintercept = 18, color="grey", linetype="solid") +
geom_vline(xintercept = 25, color="grey", linetype="solid") +
geom_point(size=2, alpha=0.7) +
geom_smooth(method="loess", level=0.95, col="#152370") +
theme_classic() +
theme(axis.title = element_text(size=24), axis.ticks.x=element_blank(), text = element_text(size = 18), axis.text.x = element_text(size=20), axis.text.y = element_text(size=22), plot.title = element_text(face="bold", hjust=0.5, size=22), legend.position="none")+
scale_shape_manual(values=c(20))+
guides(shape = FALSE, size = FALSE)+
scale_y_continuous(expand = c(0, 0), limits = c(0,0.3), breaks=seq(0,0.3,.05))+
scale_x_continuous(expand = c(0, 0), limits = c(0,26), breaks=seq(0,26,5))+
labs(y = "connectivity (Fisher z)",
x = "age (months)",
title = "")
#dev.off()
#ggsave("cingulo-operc_ant.png", path="~/Downloads", width = 4.5, height = 4)
SC_data<-subset(data_sc_1, supercluster==1 & ax=="post")
#tiff("~/Downloads/CP_2.tiff", units="in", width=4.5, height=4, res=1000)
ggplot(SC_data, aes(y=beta, x=age)) +
geom_vline(xintercept = 6, color="grey", linetype="solid") +
geom_vline(xintercept = 12, color="grey", linetype="solid") +
geom_vline(xintercept = 18, color="grey", linetype="solid") +
geom_vline(xintercept = 25, color="grey", linetype="solid") +
geom_point(size=2, alpha=0.7) +
geom_smooth(method="loess", level=0.95, col="#6577a1") +
theme_classic() +
theme(axis.title = element_text(size=24), axis.ticks.x=element_blank(), text = element_text(size = 18), axis.text.x = element_text(size=20), axis.text.y = element_text(size=22), plot.title = element_text(face="bold", hjust=0.5, size=22), legend.position="none")+
scale_shape_manual(values=c(20))+
guides(shape = FALSE, size = FALSE)+
scale_y_continuous(expand = c(0, 0), limits = c(0,0.3), breaks=seq(0,0.3,.05))+
scale_x_continuous(expand = c(0, 0), limits = c(0,26), breaks=seq(0,26,5))+
labs(y = "connectivity (Fisher z)",
x = "age (months)",
title = "")
#dev.off()
#ggsave("cingulo-operc_post.png", path="~/Downloads", width = 4.5, height = 4)
### both plots overlaid
SC_data<-subset(data_sc_1, supercluster==1)
#tiff("~/Downloads/CP_3.tiff", units="in", width=4.5, height=4, res=1000)
ggplot(SC_data, aes(y=beta, x=age, group=ax, color=ax)) +
geom_vline(xintercept = 6, color="grey", linetype="solid") +
geom_vline(xintercept = 12, color="grey", linetype="solid") +
geom_vline(xintercept = 18, color="grey", linetype="solid") +
geom_vline(xintercept = 25, color="grey", linetype="solid") +
geom_smooth(method="loess", level=0.95) +
scale_color_manual(values=c("#152370","#6577a1"))+
theme_classic() +
theme(axis.title = element_text(size=24), axis.ticks.x=element_blank(), text = element_text(size = 18), axis.text.x = element_text(size=20), axis.text.y = element_text(size=22), plot.title = element_text(face="bold", hjust=0.5, size=22), legend.position="none")+
scale_shape_manual(values=c(20))+
guides(shape = FALSE, size = FALSE)+
scale_y_continuous(expand = c(0, 0), limits = c(0,0.3), breaks=seq(0,0.3,.05))+
scale_x_continuous(expand = c(0, 0), limits = c(0,26), breaks=seq(0,26,5))+
labs(y = "connectivity (Fisher z)",
x = "age (months)",
title = "")
#dev.off()
#ggsave("cingulo-operc_antpost.png", path="~/Downloads", width = 4.5, height = 4)
SC_data<-subset(data_sc_2, supercluster==2 & ax=="ant")
#tiff("~/Downloads/FPN_1.tiff", units="in", width=4.5, height=4, res=1000)
ggplot(SC_data, aes(y=beta, x=age)) +
geom_vline(xintercept = 6, color="grey", linetype="solid") +
geom_vline(xintercept = 12, color="grey", linetype="solid") +
geom_vline(xintercept = 18, color="grey", linetype="solid") +
geom_vline(xintercept = 25, color="grey", linetype="solid") +
geom_point(size=2, alpha=0.7) +
geom_smooth(method="loess", level=0.95, col="#4bf52a") +
theme_classic() +
theme(axis.title = element_text(size=24), axis.ticks.x=element_blank(), text = element_text(size = 18), axis.text.x = element_text(size=20), axis.text.y = element_text(size=22), plot.title = element_text(face="bold", hjust=0.5, size=22), legend.position="none")+
scale_shape_manual(values=c(20))+
guides(shape = FALSE, size = FALSE)+
scale_y_continuous(expand = c(0, 0), limits = c(0,0.3), breaks=seq(0,0.3,.05))+
scale_x_continuous(expand = c(0, 0), limits = c(0,26), breaks=seq(0,26,5))+
labs(y = "connectivity (Fisher z)",
x = "age (months)",
title = "")
#dev.off()
#ggsave("DFP_ant_younger.png", path="~/Downloads", width = 4.5, height = 4)
SC_data<-subset(data_sc_2, supercluster==2 & ax=="post")
#tiff("~/Downloads/FPN_2.tiff", units="in", width=4.5, height=4, res=1000)
ggplot(SC_data, aes(y=beta, x=age)) +
geom_vline(xintercept = 6, color="grey", linetype="solid") +
geom_vline(xintercept = 12, color="grey", linetype="solid") +
geom_vline(xintercept = 18, color="grey", linetype="solid") +
geom_vline(xintercept = 25, color="grey", linetype="solid") +
geom_point(size=2, alpha=0.7) +
geom_smooth(method="loess", level=0.95, col="#a2f792") +
theme_classic() +
theme(axis.title = element_text(size=24), axis.ticks.x=element_blank(), text = element_text(size = 18), axis.text.x = element_text(size=20), axis.text.y = element_text(size=22), plot.title = element_text(face="bold", hjust=0.5, size=22), legend.position="none")+
scale_shape_manual(values=c(20))+
guides(shape = FALSE, size = FALSE)+
scale_y_continuous(expand = c(0, 0), limits = c(0,0.3), breaks=seq(0,0.3,.05))+
scale_x_continuous(expand = c(0, 0), limits = c(0,26), breaks=seq(0,26,5))+
labs(y = "connectivity (Fisher z)",
x = "age (months)",
title = "")
#dev.off()
#ggsave("DFP_post_younger.png", path="~/Downloads", width = 4.5, height = 4)
SC_data<-subset(data_sc_2, supercluster==2)
#tiff("~/Downloads/FPN_3.tiff", units="in", width=4.5, height=4, res=1000)
ggplot(SC_data, aes(y=beta, x=age, group=ax, color=ax)) +
geom_vline(xintercept = 6, color="grey", linetype="solid") +
geom_vline(xintercept = 12, color="grey", linetype="solid") +
geom_vline(xintercept = 18, color="grey", linetype="solid") +
geom_vline(xintercept = 25, color="grey", linetype="solid") +
geom_smooth(method="loess", level=0.95) +
scale_color_manual(values=c("#4bf52a","#bdf7b2"))+
theme_classic() +
theme(axis.title = element_text(size=24), axis.ticks.x=element_blank(), text = element_text(size = 18), axis.text.x = element_text(size=20), axis.text.y = element_text(size=22), plot.title = element_text(face="bold", hjust=0.5, size=22), legend.position="none")+
scale_shape_manual(values=c(20))+
guides(shape = FALSE, size = FALSE)+
scale_y_continuous(expand = c(0, 0), limits = c(0,0.3), breaks=seq(0,0.3,.05))+
scale_x_continuous(expand = c(0, 0), limits = c(0,26), breaks=seq(0,26,5))+
labs(y = "connectivity (Fisher z)",
x = "age (months)",
title = "")
#dev.off()
#ggsave("DFP_antpost.png", path="~/Downloads", width = 4.5, height = 4)
SC_data<-subset(data_sc_3, supercluster==3 & ax=="ant")
#tiff("~/Downloads/STS_1.tiff", units="in", width=4.5, height=4, res=1000)
ggplot(SC_data, aes(y=beta, x=age)) +
geom_vline(xintercept = 6, color="grey", linetype="solid") +
geom_vline(xintercept = 12, color="grey", linetype="solid") +
geom_vline(xintercept = 18, color="grey", linetype="solid") +
geom_vline(xintercept = 25, color="grey", linetype="solid") +
geom_point(size=2, alpha=0.7) +
geom_smooth(method="loess", level=0.95, col="#f2f207") +
theme_classic() +
theme(axis.title = element_text(size=24), axis.ticks.x=element_blank(), text = element_text(size = 18), axis.text.x = element_text(size=20), axis.text.y = element_text(size=22), plot.title = element_text(face="bold", hjust=0.5, size=22), legend.position="none")+
scale_shape_manual(values=c(20))+
guides(shape = FALSE, size = FALSE)+
scale_y_continuous(expand = c(0, 0), limits = c(0,0.3), breaks=seq(0,0.3,.05))+
scale_x_continuous(expand = c(0, 0), limits = c(0,26), breaks=seq(0,26,5))+
labs(y = "connectivity (Fisher z)",
x = "age (months)",
title = "")
#dev.off()
#ggsave("mPFC_ant_younger.png", path="~/Downloads", width = 4.5, height = 4)
SC_data<-subset(data_sc_3, supercluster==3 & ax=="post")
#tiff("~/Downloads/STS_2.tiff", units="in", width=4.5, height=4, res=1000)
ggplot(SC_data, aes(y=beta, x=age)) +
geom_vline(xintercept = 6, color="grey", linetype="solid") +
geom_vline(xintercept = 12, color="grey", linetype="solid") +
geom_vline(xintercept = 18, color="grey", linetype="solid") +
geom_vline(xintercept = 25, color="grey", linetype="solid") +
geom_point(size=2, alpha=0.7) +
geom_smooth(method="loess", level=0.95, col="#fafab4") +
theme_classic() +
theme(axis.title = element_text(size=24), axis.ticks.x=element_blank(), text = element_text(size = 18), axis.text.x = element_text(size=20), axis.text.y = element_text(size=22), plot.title = element_text(face="bold", hjust=0.5, size=22), legend.position="none")+
scale_shape_manual(values=c(20))+
guides(shape = FALSE, size = FALSE)+
scale_y_continuous(expand = c(0, 0), limits = c(0,0.3), breaks=seq(0,0.3,.05))+
scale_x_continuous(expand = c(0, 0), limits = c(0,26), breaks=seq(0,26,5))+
labs(y = "connectivity (Fisher z)",
x = "age (months)",
title = "")
#dev.off()
#ggsave("mPFC_post_younger.png", path="~/Downloads", width = 4.5, height = 4)
SC_data<-subset(data_sc_3, supercluster==3)
#tiff("~/Downloads/STS_3.tiff", units="in", width=4.5, height=4, res=1000)
ggplot(SC_data, aes(y=beta, x=age, group=ax, color=ax)) +
geom_vline(xintercept = 6, color="grey", linetype="solid") +
geom_vline(xintercept = 12, color="grey", linetype="solid") +
geom_vline(xintercept = 18, color="grey", linetype="solid") +
geom_vline(xintercept = 25, color="grey", linetype="solid") +
geom_smooth(method="loess", level=0.95) +
scale_color_manual(values=c("#f2f207","#fafab4"))+
theme_classic() +
theme(axis.title = element_text(size=24), axis.ticks.x=element_blank(), text = element_text(size = 18), axis.text.x = element_text(size=20), axis.text.y = element_text(size=22), plot.title = element_text(face="bold", hjust=0.5, size=22), legend.position="none")+
scale_shape_manual(values=c(20))+
guides(shape = FALSE, size = FALSE)+
scale_y_continuous(expand = c(0, 0), limits = c(0,0.3), breaks=seq(0,0.3,.05))+
scale_x_continuous(expand = c(0, 0), limits = c(0,26), breaks=seq(0,26,5))+
labs(y = "connectivity (Fisher z)",
x = "age (months)",
title = "")
#dev.off()
#ggsave("mPFC_antpost.png", path="~/Downloads", width = 4.5, height = 4)
SC_data<-subset(data_sc_5, supercluster==5 & ax=="ant")
#tiff("~/Downloads/parietal_1.tiff", units="in", width=4.5, height=4, res=1000)
ggplot(SC_data, aes(y=beta, x=age)) +
geom_vline(xintercept = 6, color="grey", linetype="solid") +
geom_vline(xintercept = 12, color="grey", linetype="solid") +
geom_vline(xintercept = 18, color="grey", linetype="solid") +
geom_vline(xintercept = 25, color="grey", linetype="solid") +
geom_point(size=2, alpha=0.7) +
geom_smooth(method="loess", level=0.95, col="#fa750f") +
theme_classic() +
theme(axis.title = element_text(size=24), axis.ticks.x=element_blank(), text = element_text(size = 18), axis.text.x = element_text(size=20), axis.text.y = element_text(size=22), plot.title = element_text(face="bold", hjust=0.5, size=22), legend.position="none")+
scale_shape_manual(values=c(20))+
guides(shape = FALSE, size = FALSE)+
scale_y_continuous(expand = c(0, 0), limits = c(0,0.3), breaks=seq(0,0.3,.05))+
scale_x_continuous(expand = c(0, 0), limits = c(0,26), breaks=seq(0,26,5))+
labs(y = "connectivity (Fisher z)",
x = "age (months)",
title = "")
#dev.off()
#ggsave("parietal_ant_younger.png", path="~/Downloads", width = 4.5, height = 4)
SC_data<-subset(data_sc_5, supercluster==5 & ax=="post")
#tiff("~/Downloads/parietal_2.tiff", units="in", width=4.5, height=4, res=1000)
ggplot(SC_data, aes(y=beta, x=age)) +
geom_vline(xintercept = 6, color="grey", linetype="solid") +
geom_vline(xintercept = 12, color="grey", linetype="solid") +
geom_vline(xintercept = 18, color="grey", linetype="solid") +
geom_vline(xintercept = 25, color="grey", linetype="solid") +
geom_point(size=2, alpha=0.7) +
geom_smooth(method="loess", level=0.95, col="#faa666") +
theme_classic() +
theme(axis.title = element_text(size=24), axis.ticks.x=element_blank(), text = element_text(size = 18), axis.text.x = element_text(size=20), axis.text.y = element_text(size=22), plot.title = element_text(face="bold", hjust=0.5, size=22), legend.position="none")+
scale_shape_manual(values=c(20))+
guides(shape = FALSE, size = FALSE)+
scale_y_continuous(expand = c(0, 0), limits = c(0,0.3), breaks=seq(0,0.3,.05))+
scale_x_continuous(expand = c(0, 0), limits = c(0,26), breaks=seq(0,26,5))+
labs(y = "connectivity (Fisher z)",
x = "age (months)",
title = "")
#dev.off()
#ggsave("parietal_post_younger.png", path="~/Downloads", width = 4.5, height = 4)
SC_data<-subset(data_sc_5, supercluster==5)
#tiff("~/Downloads/parietal_3.tiff", units="in", width=4.5, height=4, res=1000)
ggplot(SC_data, aes(y=beta, x=age, group=ax, color=ax)) +
geom_vline(xintercept = 6, color="grey", linetype="solid") +
geom_vline(xintercept = 12, color="grey", linetype="solid") +
geom_vline(xintercept = 18, color="grey", linetype="solid") +
geom_vline(xintercept = 25, color="grey", linetype="solid") +
geom_smooth(method="loess", level=0.95) +
scale_color_manual(values=c("#fa750f","#faa666"))+
theme_classic() +
theme(axis.title = element_text(size=24), axis.ticks.x=element_blank(), text = element_text(size = 18), axis.text.x = element_text(size=20), axis.text.y = element_text(size=22), plot.title = element_text(face="bold", hjust=0.5, size=22), legend.position="none")+
scale_shape_manual(values=c(20))+
guides(shape = FALSE, size = FALSE)+
scale_y_continuous(expand = c(0, 0), limits = c(0,0.3), breaks=seq(0,0.3,.05))+
scale_x_continuous(expand = c(0, 0), limits = c(0,26), breaks=seq(0,26,5))+
labs(y = "connectivity (Fisher z)",
x = "age (months)",
title = "")
#dev.off()
#ggsave("parietal_antpost.png", path="~/Downloads", width = 4.5, height = 4)
SC_data<-subset(data_sc_6, supercluster==6 & ax=="ant")
#tiff("~/Downloads/EC_1.tiff", units="in", width=4.5, height=4, res=1000)
ggplot(SC_data, aes(y=beta, x=age)) +
geom_vline(xintercept = 6, color="grey", linetype="solid") +
geom_vline(xintercept = 12, color="grey", linetype="solid") +
geom_vline(xintercept = 18, color="grey", linetype="solid") +
geom_vline(xintercept = 25, color="grey", linetype="solid") +
geom_point(size=2, alpha=0.7) +
geom_smooth(method="loess", level=0.95, col="#28fcfc") +
theme_classic() +
theme(axis.title = element_text(size=24), axis.ticks.x=element_blank(), text = element_text(size = 18), axis.text.x = element_text(size=20), axis.text.y = element_text(size=22), plot.title = element_text(face="bold", hjust=0.5, size=22), legend.position="none")+
scale_shape_manual(values=c(20))+
guides(shape = FALSE, size = FALSE)+
scale_y_continuous(expand = c(0, 0), limits = c(0,0.3), breaks=seq(0,0.3,.05))+
scale_x_continuous(expand = c(0, 0), limits = c(0,26), breaks=seq(0,26,5))+
labs(y = "connectivity (Fisher z)",
x = "age (months)",
title = "")
#dev.off()
#ggsave("entorhinal_ant_younger.png", path="~/Downloads", width = 4.5, height = 4)
SC_data<-subset(data_sc_6, supercluster==6 & ax=="post")
#tiff("~/Downloads/EC_2.tiff", units="in", width=4.5, height=4, res=1000)
ggplot(SC_data, aes(y=beta, x=age)) +
geom_vline(xintercept = 6, color="grey", linetype="solid") +
geom_vline(xintercept = 12, color="grey", linetype="solid") +
geom_vline(xintercept = 18, color="grey", linetype="solid") +
geom_vline(xintercept = 25, color="grey", linetype="solid") +
geom_point(size=2, alpha=0.7) +
geom_smooth(method="loess", level=0.95, col="#bfffff") +
theme_classic() +
theme(axis.title = element_text(size=24), axis.ticks.x=element_blank(), text = element_text(size = 18), axis.text.x = element_text(size=20), axis.text.y = element_text(size=22), plot.title = element_text(face="bold", hjust=0.5, size=22), legend.position="none")+
scale_shape_manual(values=c(20))+
guides(shape = FALSE, size = FALSE)+
scale_y_continuous(expand = c(0, 0), limits = c(0,0.3), breaks=seq(0,0.3,.05))+
scale_x_continuous(expand = c(0, 0), limits = c(0,26), breaks=seq(0,26,5))+
labs(y = "connectivity (Fisher z)",
x = "age (months)",
title = "")
#dev.off()
#ggsave("entorhinal_post_younger.png", path="~/Downloads", width = 4.5, height = 4)
SC_data<-subset(data_sc_6, supercluster==6)
#tiff("~/Downloads/EC_3.tiff", units="in", width=4.5, height=4, res=1000)
ggplot(SC_data, aes(y=beta, x=age, group=ax, color=ax)) +
geom_vline(xintercept = 6, color="grey", linetype="solid") +
geom_vline(xintercept = 12, color="grey", linetype="solid") +
geom_vline(xintercept = 18, color="grey", linetype="solid") +
geom_vline(xintercept = 25, color="grey", linetype="solid") +
geom_smooth(method="loess", level=0.95) +
scale_color_manual(values=c("#28fcfc","#bfffff"))+
theme_classic() +
theme(axis.title = element_text(size=24), axis.ticks.x=element_blank(), text = element_text(size = 18), axis.text.x = element_text(size=20), axis.text.y = element_text(size=22), plot.title = element_text(face="bold", hjust=0.5, size=22), legend.position="none")+
scale_shape_manual(values=c(20))+
guides(shape = FALSE, size = FALSE)+
scale_y_continuous(expand = c(0, 0), limits = c(0,0.3), breaks=seq(0,0.3,.05))+
scale_x_continuous(expand = c(0, 0), limits = c(0,26), breaks=seq(0,26,5))+
labs(y = "connectivity (Fisher z)",
x = "age (months)",
title = "")
#dev.off()
#ggsave("entorhinal_antpost.png", path="~/Downloads", width = 4.5, height = 4)
As the parahippocampal supercluster was the only one that did not show 100% clustering consistency, we cannot do selective averaging here. So just plot the raw data to see what it looks like. This is not corrected for subject-level bias as the other superclusters are.
## define the parahippocampal supercluster
data_SC<-subset(data_ag_hem, supercluster=="4")
## next collapse over clusters
data_SC<-aggregate(data_SC[, "beta"], by=(list(data_SC$subj, data_SC$ax, data_SC$bin)), mean)
colnames(data_SC)<-c("subj","ax","bin","beta")
## make a factor
data_SC$bin<-as.factor(data_SC$bin)
## model it
model<-lme(beta ~ ax*bin, random= ~1|subj, data=data_SC, na.action=na.omit, weights=varIdent(form=~1|bin), control = lmeControl(msMaxIter=1000, msMaxEval=1000))
## model diagnostics
# hist(residuals(model))
# plot(model)
## stats
anova(model)
## posthocs
emmeans(model, list(pairwise ~ ax|bin), adjust = "fdr")
$`emmeans of ax | bin`
bin = 1:
ax emmean SE df lower.CL upper.CL
ant 0.129 0.00695 211 0.115 0.143
post 0.222 0.00695 208 0.209 0.236
bin = 2:
ax emmean SE df lower.CL upper.CL
ant 0.168 0.00592 208 0.156 0.179
post 0.207 0.00592 208 0.195 0.218
bin = 3:
ax emmean SE df lower.CL upper.CL
ant 0.169 0.00580 208 0.158 0.180
post 0.204 0.00580 208 0.193 0.216
bin = 4:
ax emmean SE df lower.CL upper.CL
ant 0.145 0.00654 208 0.133 0.158
post 0.177 0.00654 208 0.164 0.190
Degrees-of-freedom method: containment
Confidence level used: 0.95
$`pairwise differences of ax | bin`
bin = 1:
2 estimate SE df t.ratio p.value
ant - post -0.0933 0.00583 208 -16.000 <.0001
bin = 2:
2 estimate SE df t.ratio p.value
ant - post -0.0392 0.00411 208 -9.535 <.0001
bin = 3:
2 estimate SE df t.ratio p.value
ant - post -0.0351 0.00362 208 -9.678 <.0001
bin = 4:
2 estimate SE df t.ratio p.value
ant - post -0.0318 0.00397 208 -8.013 <.0001
Degrees-of-freedom method: containment
emmeans(model, list(pairwise ~ bin|ax), adjust = "fdr")
$`emmeans of bin | ax`
ax = ant:
bin emmean SE df lower.CL upper.CL
1 0.129 0.00695 211 0.115 0.143
2 0.168 0.00592 208 0.156 0.179
3 0.169 0.00580 208 0.158 0.180
4 0.145 0.00654 208 0.133 0.158
ax = post:
bin emmean SE df lower.CL upper.CL
1 0.222 0.00695 208 0.209 0.236
2 0.207 0.00592 208 0.195 0.218
3 0.204 0.00580 208 0.193 0.216
4 0.177 0.00654 208 0.164 0.190
Degrees-of-freedom method: containment
Confidence level used: 0.95
$`pairwise differences of bin | ax`
ax = ant:
2 estimate SE df t.ratio p.value
bin1 - bin2 -0.03837 0.00913 208 -4.203 0.0001
bin1 - bin3 -0.03984 0.00905 208 -4.401 0.0001
bin1 - bin4 -0.01622 0.00954 208 -1.700 0.1087
bin2 - bin3 -0.00146 0.00828 208 -0.177 0.8601
bin2 - bin4 0.02215 0.00882 208 2.512 0.0191
bin3 - bin4 0.02361 0.00873 208 2.703 0.0149
ax = post:
2 estimate SE df t.ratio p.value
bin1 - bin2 0.01569 0.00913 208 1.719 0.1045
bin1 - bin3 0.01837 0.00905 208 2.029 0.0656
bin1 - bin4 0.04525 0.00954 208 4.742 <.0001
bin2 - bin3 0.00267 0.00828 208 0.323 0.7472
bin2 - bin4 0.02955 0.00882 208 3.352 0.0029
bin3 - bin4 0.02688 0.00873 208 3.077 0.0047
Degrees-of-freedom method: containment
P value adjustment: fdr method for 6 tests
## summarize for plotting
stats<-summarySE(data=data_SC, measurevar = "beta", groupvars = c("bin","ax"), na.rm=TRUE)
stats$upper_se<-stats$beta + stats$se
stats$lower_se<-stats$beta - stats$se
## plot
#tiff("~/Downloads/test.tiff", units="in", width=4.5, height=4, res=1000)
ggplot(stats, aes(x=bin, y=beta, group=ax)) +
geom_line(size=1, aes(linetype="solid", color=ax)) +
geom_errorbar(width = .0, aes(ymin=lower_se, ymax=upper_se)) +
theme_classic() +
theme(axis.title = element_text(size=24), axis.ticks.x=element_blank(), axis.text=element_text(size=18),
plot.title = element_text(size=26, face="bold", hjust=0.5), legend.position="bottom") +
theme(text = element_text(size = 18), axis.text.x = element_text(size=20), axis.text.y = element_text(size=22), plot.title = element_text(face="bold", hjust=0.5, size=22), legend.position="none")+
scale_x_discrete(limits=c("1", "2", "3", "4"), labels=c("0-6", "7-12", "13-18","19-25"))+
scale_y_continuous(expand = c(0, 0), limits = c(0,0.3), breaks=seq(0,0.3,.05))+
scale_color_manual(values=c('#f5200c','#fc9288')) +
geom_hline(yintercept=0)+
# ggtitle("parahippocampal
# supercluster")+
labs(x="age (months)", y="connectivity", color="hippocampal axis")
#dev.off()
#ggsave("parahippocampal.png", path="~/Downloads", width = 4.5, height = 4)
## subset
data_SC<-subset(data_ag_hem, supercluster=="4" & ax == "ant") ## it's made up of clusters 5 and 10
## next collapse over clusters
data_SC<-aggregate(data_SC[, "beta"], by=(list(data_SC$subj, data_SC$ax, data_SC$age, data_SC$supercluster)), mean)
colnames(data_SC)<-c("subj","ax", "age", "supercluster","beta")
## plot
#tiff("~/Downloads/PH_1.tiff", units="in", width=4.5, height=4, res=1000)
ggplot(data_SC, aes(y=beta, x=age)) +
geom_vline(xintercept = 6, color="grey", linetype="solid") +
geom_vline(xintercept = 12, color="grey", linetype="solid") +
geom_vline(xintercept = 18, color="grey", linetype="solid") +
geom_vline(xintercept = 25, color="grey", linetype="solid") +
geom_point(size=2, alpha=0.7) +
geom_smooth(method="loess", level=0.95, col="#f5200c") +
theme_classic() +
theme(axis.title = element_text(size=24), axis.ticks.x=element_blank(), text = element_text(size = 18), axis.text.x = element_text(size=20), axis.text.y = element_text(size=22), plot.title = element_text(face="bold", hjust=0.5, size=22), legend.position="none")+
scale_shape_manual(values=c(20))+
guides(shape = FALSE, size = FALSE)+
scale_y_continuous(expand = c(0, 0), limits = c(0,0.3), breaks=seq(0,0.3,.05))+
scale_x_continuous(expand = c(0, 0), limits = c(0,26), breaks=seq(0,26,5))+
labs(y = "connectivity (Fisher z)",
x = "age (months)",
title = "")
##dev.off()
#ggsave("parahip_ant_younger.png", path="~/Downloads", width = 4.5, height = 4)
data_SC<-subset(data_ag_hem, supercluster==4 & ax=="post")
#tiff("~/Downloads/PH_2.tiff", units="in", width=4.5, height=4, res=1000)
ggplot(data_SC, aes(y=beta, x=age)) +
geom_vline(xintercept = 6, color="grey", linetype="solid") +
geom_vline(xintercept = 12, color="grey", linetype="solid") +
geom_vline(xintercept = 18, color="grey", linetype="solid") +
geom_vline(xintercept = 25, color="grey", linetype="solid") +
geom_point(size=2, alpha=0.7) +
geom_smooth(method="loess", level=0.95, col="#fc9288") +
theme_classic() +
theme(axis.title = element_text(size=24), axis.ticks.x=element_blank(), text = element_text(size = 18), axis.text.x = element_text(size=20), axis.text.y = element_text(size=22), plot.title = element_text(face="bold", hjust=0.5, size=22), legend.position="none")+
scale_shape_manual(values=c(20))+
guides(shape = FALSE, size = FALSE)+
scale_y_continuous(expand = c(0, 0), limits = c(0,0.3), breaks=seq(0,0.3,.05))+
scale_x_continuous(expand = c(0, 0), limits = c(0,26), breaks=seq(0,26,5))+
labs(y = "connectivity (Fisher z)",
x = "age (months)",
title = "")
#dev.off()
#ggsave("parahip_post_younger.png", path="~/Downloads", width = 4.5, height = 4)
### both plots overlaid
data_SC<-subset(data_ag_hem, supercluster==4)
#tiff("~/Downloads/PH_3.tiff", units="in", width=4.5, height=4, res=1000)
ggplot(data_SC, aes(y=beta, x=age, group=ax, color=ax)) +
geom_vline(xintercept = 6, color="grey", linetype="solid") +
geom_vline(xintercept = 12, color="grey", linetype="solid") +
geom_vline(xintercept = 18, color="grey", linetype="solid") +
geom_vline(xintercept = 25, color="grey", linetype="solid") +
geom_smooth(method="loess", level=0.95) +
scale_color_manual(values=c("#f5200c","#fc9288"))+
theme_classic() +
theme(axis.title = element_text(size=24), axis.ticks.x=element_blank(), text = element_text(size = 18), axis.text.x = element_text(size=20), axis.text.y = element_text(size=22), plot.title = element_text(face="bold", hjust=0.5, size=22), legend.position="none")+
scale_shape_manual(values=c(20))+
guides(shape = FALSE, size = FALSE)+
scale_y_continuous(expand = c(0, 0), limits = c(0,0.3), breaks=seq(0,0.3,.05))+
scale_x_continuous(expand = c(0, 0), limits = c(0,26), breaks=seq(0,26,5))+
labs(y = "connectivity (Fisher z)",
x = "age (months)",
title = "")
#dev.off()
#ggsave("parahip_antpost.png", path="~/Downloads", width = 4.5, height = 4)
Comparing connectivity of the anterior and posterior hippocampus to each supercluster, between the oldest age bin within the infantile amnesia window (~1.5-2 years old) and children just outside of it (3-6 years)
## read in the data, CHANGE PATHS HERE
dir<-"/Users/samanthaaudrain/Documents/Science_Projects/BabyHippos/R_BabyHippos"
dir<-"/Users/audrainsp/Library/CloudStorage/OneDrive-NationalInstitutesofHealth/BabyHippos/R_BabyHippos"
data2 <- read.table(paste(dir, "/supercluster_antpost_betas_by_age.txt", sep=""), header = TRUE)
data2<- subset(data2, age > 35)
data2$bin<-"5"
## change roi label from ant_hippo to ant and post to be consistent with other data
data2<-mutate(data2, ax= case_when(roi=="ant_hippo" ~ "ant",
roi=="post_hippo"~ "post"))
## collapse over hem for each subj
data2_agg<-aggregate(data2[, "beta"], by=(list(data2$subj, data2$bin, data2$age, data2$ax, data2$supercluster)), mean)
colnames(data2_agg)<-c("subj","bin","age","ax","supercluster","beta")
## reorder
data2_agg<-as.data.frame(cbind(data2_agg$subj, data2_agg$bin, data2_agg$age, data2_agg$ax, data2_agg$beta, data2_agg$supercluster))
colnames(data2_agg)<-c("subj","bin","age","ax","beta","supercluster")
## remove parahip cluster which doesnt have LOO data to compare to
data2_agg_nosc4<-subset(data2_agg, supercluster != 4)
## combine original LOO dataset with the data from the older kids
data_LOO<-rbind(data_SC_all, data2_agg_nosc4)
## make age numeric and subset the data to be just the oldest age bin within infantile amnesia window + the older kids
data_LOO$age<-as.numeric(data_LOO$age)
data_LOO<-subset(data_LOO, ((age>18 & age<26) | age>35))
## specify other variables as factors or numeric
data_LOO$supercluster<-as.factor(data_LOO$supercluster)
data_LOO$beta <- as.numeric(data_LOO$beta)
## create labels for inside and outside infantile amnesia window
data_LOO<-data_LOO %>%
mutate(amnesia=case_when(age<26 ~ "inside",
age>35 ~ "outside"))
## model
model<- lme((beta) ~ ax*amnesia*supercluster, random=~supercluster|subj, weights=varIdent(form=~1|supercluster), data=data_LOO, control = lmeControl(msMaxIter=1000, msMaxEval=1000), na.action=na.omit)
## modeled heteroskedasticity
## model diagnostics
hist(residuals(model))
plot(model)
## stats
anova(model)
## post-hocs
emmeans(model, list(pairwise ~ amnesia | supercluster), adjust = "fdr")
NOTE: Results may be misleading due to involvement in interactions
$`emmeans of amnesia | supercluster`
supercluster = 1:
amnesia emmean SE df lower.CL upper.CL
inside 0.0529 0.00526 73 0.0424 0.0633
outside 0.1298 0.00655 72 0.1168 0.1429
supercluster = 2:
amnesia emmean SE df lower.CL upper.CL
inside 0.0463 0.00624 73 0.0338 0.0587
outside 0.0867 0.00777 72 0.0712 0.1022
supercluster = 3:
amnesia emmean SE df lower.CL upper.CL
inside 0.0857 0.00685 73 0.0720 0.0993
outside 0.1064 0.00853 72 0.0894 0.1234
supercluster = 5:
amnesia emmean SE df lower.CL upper.CL
inside 0.1000 0.00635 73 0.0873 0.1127
outside 0.1607 0.00792 72 0.1449 0.1764
supercluster = 6:
amnesia emmean SE df lower.CL upper.CL
inside 0.1661 0.00732 73 0.1515 0.1807
outside 0.1606 0.00912 72 0.1424 0.1788
Results are averaged over the levels of: ax
Degrees-of-freedom method: containment
Results are given on the ( (not the response) scale.
Confidence level used: 0.95
$`pairwise differences of amnesia | supercluster`
supercluster = 1:
2 estimate SE df t.ratio p.value
inside - outside -0.07697 0.00840 72 -9.162 <.0001
supercluster = 2:
2 estimate SE df t.ratio p.value
inside - outside -0.04046 0.00997 72 -4.060 0.0001
supercluster = 3:
2 estimate SE df t.ratio p.value
inside - outside -0.02074 0.01090 72 -1.896 0.0620
supercluster = 5:
2 estimate SE df t.ratio p.value
inside - outside -0.06067 0.01020 72 -5.977 <.0001
supercluster = 6:
2 estimate SE df t.ratio p.value
inside - outside 0.00548 0.01170 72 0.469 0.6407
Results are averaged over the levels of: ax
Note: contrasts are still on the ( scale. Consider using
regrid() if you want contrasts of back-transformed estimates.
Degrees-of-freedom method: containment
emmeans(model, list(pairwise ~ ax*supercluster), adjust = "fdr")
NOTE: Results may be misleading due to involvement in interactions
$`emmeans of ax, supercluster`
ax supercluster emmean SE df lower.CL upper.CL
ant 1 0.0832 0.00438 72 0.0745 0.0920
post 1 0.0994 0.00438 72 0.0907 0.1082
ant 2 0.0543 0.00514 72 0.0441 0.0646
post 2 0.0787 0.00514 72 0.0684 0.0889
ant 3 0.0987 0.00569 72 0.0874 0.1101
post 3 0.0933 0.00569 72 0.0820 0.1047
ant 5 0.1236 0.00547 72 0.1127 0.1345
post 5 0.1371 0.00547 72 0.1262 0.1480
ant 6 0.1891 0.00644 72 0.1763 0.2020
post 6 0.1376 0.00644 72 0.1247 0.1504
Results are averaged over the levels of: amnesia
Degrees-of-freedom method: containment
Results are given on the ( (not the response) scale.
Confidence level used: 0.95
$`pairwise differences of ax, supercluster`
1 estimate SE df t.ratio p.value
ant supercluster1 - post supercluster1 -0.016170 0.00249 648 -6.485 <.0001
ant supercluster1 - ant supercluster2 0.028917 0.00335 648 8.635 <.0001
ant supercluster1 - post supercluster2 0.004599 0.00335 648 1.373 0.1823
ant supercluster1 - ant supercluster3 -0.015495 0.00454 648 -3.412 0.0009
ant supercluster1 - post supercluster3 -0.010096 0.00454 648 -2.223 0.0298
ant supercluster1 - ant supercluster5 -0.040308 0.00467 648 -8.625 <.0001
ant supercluster1 - post supercluster5 -0.053852 0.00467 648 -11.523 <.0001
ant supercluster1 - ant supercluster6 -0.105876 0.00569 648 -18.602 <.0001
ant supercluster1 - post supercluster6 -0.054306 0.00569 648 -9.542 <.0001
post supercluster1 - ant supercluster2 0.045087 0.00335 648 13.464 <.0001
post supercluster1 - post supercluster2 0.020769 0.00335 648 6.202 <.0001
post supercluster1 - ant supercluster3 0.000675 0.00454 648 0.149 0.9019
post supercluster1 - post supercluster3 0.006074 0.00454 648 1.338 0.1899
post supercluster1 - ant supercluster5 -0.024138 0.00467 648 -5.165 <.0001
post supercluster1 - post supercluster5 -0.037682 0.00467 648 -8.063 <.0001
post supercluster1 - ant supercluster6 -0.089705 0.00569 648 -15.761 <.0001
post supercluster1 - post supercluster6 -0.038136 0.00569 648 -6.700 <.0001
ant supercluster2 - post supercluster2 -0.024318 0.00253 648 -9.610 <.0001
ant supercluster2 - ant supercluster3 -0.044412 0.00442 648 -10.047 <.0001
ant supercluster2 - post supercluster3 -0.039013 0.00442 648 -8.826 <.0001
ant supercluster2 - ant supercluster5 -0.069226 0.00485 648 -14.282 <.0001
ant supercluster2 - post supercluster5 -0.082769 0.00485 648 -17.076 <.0001
ant supercluster2 - ant supercluster6 -0.134793 0.00615 648 -21.926 <.0001
ant supercluster2 - post supercluster6 -0.083224 0.00615 648 -13.537 <.0001
post supercluster2 - ant supercluster3 -0.020094 0.00442 648 -4.546 <.0001
post supercluster2 - post supercluster3 -0.014695 0.00442 648 -3.324 0.0011
post supercluster2 - ant supercluster5 -0.044907 0.00485 648 -9.265 <.0001
post supercluster2 - post supercluster5 -0.058451 0.00485 648 -12.059 <.0001
post supercluster2 - ant supercluster6 -0.110475 0.00615 648 -17.970 <.0001
post supercluster2 - post supercluster6 -0.058905 0.00615 648 -9.582 <.0001
ant supercluster3 - post supercluster3 0.005399 0.00313 648 1.723 0.0937
ant supercluster3 - ant supercluster5 -0.024813 0.00518 648 -4.786 <.0001
ant supercluster3 - post supercluster5 -0.038357 0.00518 648 -7.399 <.0001
ant supercluster3 - ant supercluster6 -0.090380 0.00630 648 -14.341 <.0001
ant supercluster3 - post supercluster6 -0.038811 0.00630 648 -6.158 <.0001
post supercluster3 - ant supercluster5 -0.030212 0.00518 648 -5.828 <.0001
post supercluster3 - post supercluster5 -0.043756 0.00518 648 -8.440 <.0001
post supercluster3 - ant supercluster6 -0.095779 0.00630 648 -15.198 <.0001
post supercluster3 - post supercluster6 -0.044210 0.00630 648 -7.015 <.0001
ant supercluster5 - post supercluster5 -0.013544 0.00407 648 -3.330 0.0011
ant supercluster5 - ant supercluster6 -0.065567 0.00616 648 -10.645 <.0001
ant supercluster5 - post supercluster6 -0.013998 0.00616 648 -2.273 0.0270
post supercluster5 - ant supercluster6 -0.052024 0.00616 648 -8.446 <.0001
post supercluster5 - post supercluster6 -0.000454 0.00616 648 -0.074 0.9412
ant supercluster6 - post supercluster6 0.051569 0.00537 648 9.601 <.0001
Results are averaged over the levels of: amnesia
Note: contrasts are still on the ( scale. Consider using
regrid() if you want contrasts of back-transformed estimates.
Degrees-of-freedom method: containment
P value adjustment: fdr method for 45 tests
The model above was the best fit and had the best diagnostics, these weren’t used in the end
# anova(model2,model)
# model<-lme((beta) ~ ax*amnesia*supercluster, random=~supercluster*amnesia|subj, weights=varIdent(form=~1|supercluster), data=data_ag_hem, control = lmeControl(msMaxIter=1000, msMaxEval=1000), na.action=na.omit) ## doesnt do better than version where just have supercluster random slope. cant have 3 way random slope.
# model<-lme(beta ~ ax*amnesia*supercluster, random=~1|subj, weights=varIdent(form=~1|amnesia*ax*supercluster), data=data_ag_hem, control = lmeControl(msMaxIter=1000, msMaxEval=1000), na.action=na.omit) ## original model, heteroskedastic and nonlinear and not normal
# model<-lme(beta ~ ax*amnesia*supercluster, random=~1|subj, weights=varIdent(form=~fitted(.)), data=data_ag_hem, control = lmeControl(msMaxIter=1000, msMaxEval=1000), na.action=na.omit) ## models variance as a function of the fitted values, didnt help
# model<-lme(sqrt(beta+1) ~ ax*amnesia*supercluster, random=~1|subj, data=data_ag_hem, control = lmeControl(msMaxIter=1000, msMaxEval=1000), na.action=na.omit) ## squrt transform to help with skew didnt help
# model<-lme(log(beta+1) ~ ax*amnesia*supercluster, random=~1|subj, data=data_ag_hem, control = lmeControl(msMaxIter=1000, msMaxEval=1000), na.action=na.omit) ### log to help with skew didn thelp
# model<-lme(beta ~ ax*amnesia*supercluster, random=~supercluster|subj, data=data_ag_hem, control = lmeControl(msMaxIter=1000, msMaxEval=1000), na.action=na.omit) ## this one helped, adding a random slope for supercluster. But still supercluster a bit heteroskedastic, try adding weights
### diagnostics checked
# ## residuals should be normally distributed
# hist(residuals(model))
# ##check for homogeneity of variance
# plot(model) #should look fairly uniform and random
# leveneTest(residuals(model) ~ data_SC$ax)
# leveneTest(residuals(model) ~ data_SC$amnesia)
# shapiro.test(resid(model))
# icc(model)
# ## A high ICC suggests that subjects account for a large proportion of variance in beta, meaning subject-level variation is important to model (random intercept for subject). ICC>.2 justifies using a random intercept, if greater than .5 most of the variance is due to subject differences and a mixed model is crucial. If ICC is moderate 0.05-.2 a random intercept model is fine, if ICC is greater than .2 should consider addign random slopes if suspect effect of fixed effects varies by subject. Want to use adjusted ICC.
# AIC(model,model2) ## lower, more neg is better
# ## check variance structure of each level of amnesia
# plot(fitted(model), resid(model),
# col=as.numeric(as.factor(data_SC$amnesia)), # Convert factor to numeric for colors
# pch=19)
# abline(h=0, lty=2)
# ## if residuals look unequal, variance structure is needed
# boxplot(resid(model) ~ amnesia, data=data_SC)
## next break down model into 2 way interaction. i.e. ax by amnesia for each supercluster
## these models were the best fits after some model comparisons
data_SC<-subset(data_LOO, supercluster=="1")
model<- lme((beta) ~ ax*amnesia, random=~amnesia|subj, weights=varIdent(form=~1|ax), data=data_SC, control = lmeControl(msMaxIter=1000, msMaxEval=1000), na.action=na.omit)
anova(model)
emmeans(model, list(pairwise ~ amnesia*ax), adjust = "fdr")
$`emmeans of amnesia, ax`
amnesia ax emmean SE df lower.CL upper.CL
inside ant 0.0458 0.00398 73 0.0379 0.0538
outside ant 0.1207 0.00839 72 0.1039 0.1374
inside post 0.0599 0.00457 72 0.0507 0.0690
outside post 0.1390 0.00885 72 0.1214 0.1566
Degrees-of-freedom method: containment
Results are given on the ( (not the response) scale.
Confidence level used: 0.95
$`pairwise differences of amnesia, ax`
1 estimate SE df t.ratio p.value
inside ant - outside ant -0.0748 0.00929 72 -8.052 <.0001
inside ant - inside post -0.0140 0.00312 72 -4.486 <.0001
inside ant - outside post -0.0931 0.00970 72 -9.600 <.0001
outside ant - inside post 0.0608 0.00956 72 6.361 <.0001
outside ant - outside post -0.0183 0.00389 72 -4.715 <.0001
inside post - outside post -0.0791 0.00996 72 -7.947 <.0001
Note: contrasts are still on the ( scale. Consider using
regrid() if you want contrasts of back-transformed estimates.
Degrees-of-freedom method: containment
P value adjustment: fdr method for 6 tests
data_SC<-subset(data_LOO, supercluster=="2")
model<- lme((beta) ~ ax*amnesia, random=~amnesia|subj, weights=varIdent(form=~1|ax), data=data_SC, control = lmeControl(msMaxIter=1000, msMaxEval=1000), na.action=na.omit)
anova(model)
emmeans(model, list(pairwise ~ amnesia*ax), adjust = "fdr")
$`emmeans of amnesia, ax`
amnesia ax emmean SE df lower.CL upper.CL
inside ant 0.0335 0.00533 73 0.0229 0.0441
outside ant 0.0751 0.00861 72 0.0580 0.0923
inside post 0.0590 0.00620 72 0.0467 0.0714
outside post 0.0983 0.00947 72 0.0794 0.1172
Degrees-of-freedom method: containment
Results are given on the ( (not the response) scale.
Confidence level used: 0.95
$`pairwise differences of amnesia, ax`
1 estimate SE df t.ratio p.value
inside ant - outside ant -0.0416 0.01010 72 -4.111 0.0002
inside ant - inside post -0.0255 0.00317 72 -8.046 <.0001
inside ant - outside post -0.0648 0.01090 72 -5.960 <.0001
outside ant - inside post 0.0161 0.01060 72 1.521 0.1326
outside ant - outside post -0.0231 0.00395 72 -5.865 <.0001
inside post - outside post -0.0393 0.01130 72 -3.470 0.0011
Note: contrasts are still on the ( scale. Consider using
regrid() if you want contrasts of back-transformed estimates.
Degrees-of-freedom method: containment
P value adjustment: fdr method for 6 tests
data_SC<-subset(data_LOO, supercluster=="3")
model<- lme((beta) ~ ax*amnesia, random=~amnesia|subj, weights=varIdent(form=~1|ax), data=data_SC, control = lmeControl(msMaxIter=1000, msMaxEval=1000), na.action=na.omit)
anova(model)
emmeans(model, list(pairwise ~ amnesia*ax), adjust = "fdr")
$`emmeans of amnesia, ax`
amnesia ax emmean SE df lower.CL upper.CL
inside ant 0.0906 0.00596 73 0.0787 0.102
outside ant 0.1069 0.00979 72 0.0874 0.126
inside post 0.0808 0.00661 72 0.0676 0.094
outside post 0.1059 0.01040 72 0.0851 0.127
Degrees-of-freedom method: containment
Results are given on the ( (not the response) scale.
Confidence level used: 0.95
$`pairwise differences of amnesia, ax`
1 estimate SE df t.ratio p.value
inside ant - outside ant -0.01635 0.01150 72 -1.427 0.2370
inside ant - inside post 0.00979 0.00392 72 2.495 0.0893
inside ant - outside post -0.01534 0.01200 72 -1.278 0.2462
outside ant - inside post 0.02614 0.01180 72 2.213 0.0902
outside ant - outside post 0.00101 0.00489 72 0.207 0.8367
inside post - outside post -0.02513 0.01230 72 -2.037 0.0906
Note: contrasts are still on the ( scale. Consider using
regrid() if you want contrasts of back-transformed estimates.
Degrees-of-freedom method: containment
P value adjustment: fdr method for 6 tests
data_SC<-subset(data_LOO, supercluster=="5")
model<- lme((beta) ~ ax*amnesia, random=~amnesia|subj, weights=varComb(varIdent(form=~1|ax), varPower(form=~fitted(.))), data=data_SC, control = lmeControl(msMaxIter=1000, msMaxEval=1000), na.action=na.omit)
anova(model)
emmeans(model, list(pairwise ~ ax), adjust = "fdr")
NOTE: Results may be misleading due to involvement in interactions
$`emmeans of ax`
ax emmean SE df lower.CL upper.CL
ant 0.123 0.00536 72 0.112 0.133
post 0.136 0.00622 72 0.124 0.149
Results are averaged over the levels of: amnesia
Degrees-of-freedom method: containment
Results are given on the ( (not the response) scale.
Confidence level used: 0.95
$`pairwise differences of ax`
1 estimate SE df t.ratio p.value
ant - post -0.0135 0.00411 72 -3.282 0.0016
Results are averaged over the levels of: amnesia
Note: contrasts are still on the ( scale. Consider using
regrid() if you want contrasts of back-transformed estimates.
Degrees-of-freedom method: containment
data_SC<-subset(data_LOO, supercluster=="6")
model<- lme((beta) ~ ax*amnesia, random=~amnesia|subj, weights=varIdent(form=~1|ax), data=data_SC, control = lmeControl(msMaxIter=1000, msMaxEval=1000), na.action=na.omit)
anova(model)
emmeans(model, list(pairwise ~ amnesia*ax), adjust = "fdr")
$`emmeans of amnesia, ax`
amnesia ax emmean SE df lower.CL upper.CL
inside ant 0.194 0.00889 73 0.176 0.211
outside ant 0.185 0.00942 72 0.166 0.203
inside post 0.139 0.00820 72 0.122 0.155
outside post 0.137 0.00839 72 0.120 0.153
Degrees-of-freedom method: containment
Results are given on the ( (not the response) scale.
Confidence level used: 0.95
$`pairwise differences of amnesia, ax`
1 estimate SE df t.ratio p.value
inside ant - outside ant 0.00892 0.01300 72 0.688 0.5921
inside ant - inside post 0.05500 0.00672 72 8.179 <.0001
inside ant - outside post 0.05705 0.01220 72 4.667 <.0001
outside ant - inside post 0.04609 0.01250 72 3.691 0.0006
outside ant - outside post 0.04813 0.00838 72 5.746 <.0001
inside post - outside post 0.00205 0.01170 72 0.175 0.8619
Note: contrasts are still on the ( scale. Consider using
regrid() if you want contrasts of back-transformed estimates.
Degrees-of-freedom method: containment
P value adjustment: fdr method for 6 tests
## define the cluster
data_SC<-subset(data_LOO, supercluster=="3")
## summarize for plotting
stats<-summarySE(data=data_SC, measurevar = "beta", groupvars = c("amnesia","ax"), na.rm=TRUE)
stats$upper_se<-stats$beta + stats$se
stats$lower_se<-stats$beta - stats$se
## plot
#tiff("~/Downloads/windows_STS.tiff", units="in", width=5, height=4.5, res=1000)
ggplot(stats, aes(x=amnesia, y=beta, fill=ax)) +
theme_classic() +
geom_bar(stat = "identity", position = "dodge", width = 0.8) +
geom_errorbar(width = .0, position = position_dodge(0.8), aes(ymin=lower_se, ymax=upper_se))+
labs(x="", y="connectivity (Fisher z)") +
theme(axis.title = element_text(size=24), axis.ticks.x=element_blank(), text = element_text(size = 18), axis.text.x = element_text(size=18), axis.text.y = element_text(size=22), plot.title = element_text(face="bold", hjust=0.5, size=22), legend.position="none")+
labs(fill="") +
scale_fill_manual(values=c('#f2f207','#fafab4'), labels = c("anterior", "posterior"))+ ### mPFC-STS
scale_x_discrete(limits=c("inside", "outside"), labels=c("inside window", "outside window"))+
scale_y_continuous(expand = c(0, 0), limits = c(0,0.3), breaks=seq(0,0.3,.05))+
geom_hline(yintercept=0)+
ggtitle("vmPFC-STS
supercluster")+
theme(plot.title = element_text(size=22, face="bold", hjust=0.5), legend.position="bottom")
#dev.off()
#ggsave("vmpfc_sts.png", path="~/Downloads", width = 4.5, height = 4.5)
## define the cluster
data_SC<-subset(data_LOO, supercluster=="1")
## summarize for plotting
stats<-summarySE(data=data_SC, measurevar = "beta", groupvars = c("amnesia","ax"), na.rm=TRUE)
stats$upper_se<-stats$beta + stats$se
stats$lower_se<-stats$beta - stats$se
## plot
#tiff("~/Downloads/windows_CO.tiff", units="in", width=5, height=4.5, res=1000)
ggplot(stats, aes(x=amnesia, y=beta, fill=ax)) +
theme_classic() +
geom_bar(stat = "identity", position = "dodge", width = 0.8) +
geom_errorbar(width = .0, position = position_dodge(0.8), aes(ymin=lower_se, ymax=upper_se))+
labs(x="", y="connectivity") +
theme(axis.title = element_text(size=24), axis.ticks.x=element_blank(), text = element_text(size = 18), axis.text.x = element_text(size=18), axis.text.y = element_text(size=22), plot.title = element_text(face="bold", hjust=0.5, size=22), legend.position="none")+
labs(fill="") +
scale_fill_manual(values=c('#152370','#6577a1'), labels = c("anterior", "posterior"))+
scale_x_discrete(limits=c("inside", "outside"), labels=c("inside window", "outside window"))+
scale_y_continuous(expand = c(0, 0), limits = c(0,0.3), breaks=seq(0,0.3,.05))+
geom_hline(yintercept=0)+
ggtitle("cingulo-opercular
supercluster")+
theme(plot.title = element_text(size=22, face="bold", hjust=0.5), legend.position="bottom")
#dev.off()
#ggsave("cingulo-operc.png", path="~/Downloads", width = 4.5, height = 4.5)
## define the cluster
data_SC<-subset(data_LOO, supercluster=="2")
## summarize for plotting
stats<-summarySE(data=data_SC, measurevar = "beta", groupvars = c("amnesia","ax"), na.rm=TRUE)
stats$upper_se<-stats$beta + stats$se
stats$lower_se<-stats$beta - stats$se
## plot
#tiff("~/Downloads/windows_dFP.tiff", units="in", width=5, height=4.5, res=1000)
ggplot(stats, aes(x=amnesia, y=beta, fill=ax)) +
theme_classic() +
geom_bar(stat = "identity", position = "dodge", width = 0.8) +
geom_errorbar(width = .0, position = position_dodge(0.8), aes(ymin=lower_se, ymax=upper_se))+
labs(x="", y="connectivity (Fisher z)") +
theme(axis.title = element_text(size=24), axis.ticks.x=element_blank(), text = element_text(size = 18), axis.text.x = element_text(size=18), axis.text.y = element_text(size=22), plot.title = element_text(face="bold", hjust=0.5, size=22), legend.position="none")+
labs(fill="") +
scale_fill_manual(values=c('#4bf52a','#a2f792'), labels = c("anterior", "posterior"))+
scale_x_discrete(limits=c("inside", "outside"), labels=c("inside window", "outside window"))+
scale_y_continuous(expand = c(0, 0), limits = c(0,0.3), breaks=seq(0,0.3,.05))+
geom_hline(yintercept=0)+
ggtitle("dorsal frontal parietal
supercluster")+
theme(plot.title = element_text(size=22, face="bold", hjust=0.5), legend.position="bottom")
#dev.off()
#ggsave("dorsalFPN.png", path="~/Downloads", width = 4.5, height = 4.5)
## define the cluster
data_SC<-subset(data_LOO, supercluster=="6")
## summarize for plotting
stats<-summarySE(data=data_SC, measurevar = "beta", groupvars = c("amnesia","ax"), na.rm=TRUE)
stats$upper_se<-stats$beta + stats$se
stats$lower_se<-stats$beta - stats$se
## plot
#tiff("~/Downloads/windows_EC.tiff", units="in", width=5, height=4.5, res=1000)
ggplot(stats, aes(x=amnesia, y=beta, fill=ax)) +
theme_classic() +
geom_bar(stat = "identity", position = "dodge", width = 0.8) +
geom_errorbar(width = .0, position = position_dodge(0.8), aes(ymin=lower_se, ymax=upper_se))+
labs(x="", y="connectivity (Fisher z)") +
theme(axis.title = element_text(size=24), axis.ticks.x=element_blank(), text = element_text(size = 18), axis.text.x = element_text(size=18), axis.text.y = element_text(size=22), plot.title = element_text(face="bold", hjust=0.5, size=22), legend.position="none")+
labs(fill="") +
scale_fill_manual(values=c('#28fcfc','#bfffff'), labels = c("anterior", "posterior"))+
scale_x_discrete(limits=c("inside", "outside"), labels=c("inside window", "outside window"))+
scale_y_continuous(expand = c(0, 0), limits = c(0,0.3), breaks=seq(0,0.3,.05))+
geom_hline(yintercept=0)+
ggtitle("entorhinal
supercluster")+
theme(plot.title = element_text(size=22, face="bold", hjust=0.5), legend.position="bottom")
#dev.off()
#ggsave("entorhinal.png", path="~/Downloads", width = 4.5, height = 4.5)
## define the cluster
data_SC<-subset(data_LOO, supercluster=="5")
## summarize for plotting
stats<-summarySE(data=data_SC, measurevar = "beta", groupvars = c("amnesia","ax"), na.rm=TRUE)
stats$upper_se<-stats$beta + stats$se
stats$lower_se<-stats$beta - stats$se
## plot
#tiff("~/Downloads/windows_parietal.tiff", units="in", width=5, height=4.5, res=1000)
ggplot(stats, aes(x=amnesia, y=beta, fill=ax)) +
theme_classic() +
geom_bar(stat = "identity", position = "dodge", width = 0.8) +
geom_errorbar(width = .0, position = position_dodge(0.8), aes(ymin=lower_se, ymax=upper_se))+
labs(x="", y="connectivity") +
theme(axis.title = element_text(size=24), axis.ticks.x=element_blank(), text = element_text(size = 18), axis.text.x = element_text(size=18), axis.text.y = element_text(size=22), plot.title = element_text(face="bold", hjust=0.5, size=22), legend.position="none")+
labs(fill="") +
scale_fill_manual(values=c('#fa750f','#faa666'), labels = c("anterior", "posterior"))+
scale_x_discrete(limits=c("inside", "outside"), labels=c("inside window", "outside window"))+
scale_y_continuous(expand = c(0, 0), limits = c(0,0.3), breaks=seq(0,0.3,.05))+
geom_hline(yintercept=0)+
ggtitle("medial parietal
supercluster")+
theme(plot.title = element_text(size=22, face="bold", hjust=0.5), legend.position="bottom")
#dev.off()
#ggsave("medialparietal.png", path="~/Downloads", width = 4.5, height = 4.5)
This cluster is not part of LOO analysis so running stats and plot separately
## read in the data, CHANGE TO YOUR PATH
dir<-"/Users/audrainsp/Library/CloudStorage/OneDrive-NationalInstitutesofHealth/BabyHippos/R_BabyHippos"
data <- read.table(paste(dir, "/supercluster_antpost_betas_by_age.txt", sep=""), header = TRUE)
## get correct age range
data<-subset(data, ((age>18 & age<26) | age>35))
## collapse over hem
data_ag_hem<-aggregate(data[, "beta"], by=(list(data$subj, data$roi, data$age, data$supercluster)), mean)
colnames(data_ag_hem)<-c("subj","ax","age","supercluster","beta")
## define supercluster
data_SC<-subset(data_ag_hem, supercluster=="4")
## create lable for inside and outside infantile amnesia window
data_SC<-data_SC %>%
mutate(amnesia=case_when(age<26 ~ "inside",
age>35 ~ "outside"))
## model
model<- lme((beta) ~ ax*amnesia, random=~amnesia|subj, weights=varComb(varIdent(form=~1|amnesia), varPower(form=~fitted(.))), data=data_SC, control = lmeControl(msMaxIter=1000, msMaxEval=1000), na.action=na.omit)
## diagnostics
# hist(residuals(model))
# ##check for homogeneity of variance
# plot(model) #should look fairly uniform and random
# qqnorm(resid(model))
# qqline(resid(model))
# icc(model)
## stats
anova(model)
## posthocs
emmeans(model, list(pairwise ~ amnesia*ax), adjust = "fdr")
$`emmeans of amnesia, ax`
amnesia ax emmean SE df lower.CL upper.CL
inside ant_hippo 0.143 0.00552 73 0.132 0.154
outside ant_hippo 0.167 0.00944 72 0.148 0.186
inside post_hippo 0.175 0.00578 72 0.163 0.186
outside post_hippo 0.199 0.00961 72 0.180 0.219
Degrees-of-freedom method: containment
Results are given on the ( (not the response) scale.
Confidence level used: 0.95
$`pairwise differences of amnesia, ax`
1 estimate SE df t.ratio p.value
inside ant_hippo - outside ant_hippo -0.02355 0.01090 72 -2.154 0.0415
inside ant_hippo - inside post_hippo -0.03130 0.00352 72 -8.883 <.0001
inside ant_hippo - outside post_hippo -0.05616 0.01110 72 -5.064 <.0001
outside ant_hippo - inside post_hippo -0.00775 0.01110 72 -0.700 0.4862
outside ant_hippo - outside post_hippo -0.03260 0.00397 72 -8.215 <.0001
inside post_hippo - outside post_hippo -0.02486 0.01120 72 -2.216 0.0415
Note: contrasts are still on the ( scale. Consider using
regrid() if you want contrasts of back-transformed estimates.
Degrees-of-freedom method: containment
P value adjustment: fdr method for 6 tests
emmeans(model, list(pairwise ~ amnesia), adjust = "fdr")
NOTE: Results may be misleading due to involvement in interactions
$`emmeans of amnesia`
amnesia emmean SE df lower.CL upper.CL
inside 0.159 0.00537 72 0.148 0.170
outside 0.183 0.00932 72 0.165 0.202
Results are averaged over the levels of: ax
Degrees-of-freedom method: containment
Results are given on the ( (not the response) scale.
Confidence level used: 0.95
$`pairwise differences of amnesia`
1 estimate SE df t.ratio p.value
inside - outside -0.0242 0.0108 72 -2.251 0.0275
Results are averaged over the levels of: ax
Note: contrasts are still on the ( scale. Consider using
regrid() if you want contrasts of back-transformed estimates.
Degrees-of-freedom method: containment
## format for plotting
stats<-summarySE(data=data_SC, measurevar = "beta", groupvars = c("amnesia","ax"), na.rm=TRUE)
stats$upper_se<-stats$beta + stats$se
stats$lower_se<-stats$beta - stats$se
## plot
#tiff("~/Downloads/windows_PH.tiff", units="in", width=5, height=4.5, res=1000)
ggplot(stats, aes(x=amnesia, y=beta, fill=ax)) +
theme_classic() +
geom_bar(stat = "identity", position = "dodge", width = 0.8) +
geom_errorbar(width = .0, position = position_dodge(0.8), aes(ymin=lower_se, ymax=upper_se))+
labs(x="", y="connectivity (Fisher z)") +
theme(axis.title = element_text(size=24), axis.ticks.x=element_blank(), text = element_text(size = 18), axis.text.x = element_text(size=18), axis.text.y = element_text(size=22), plot.title = element_text(face="bold", hjust=0.5, size=22), legend.position="none")+
labs(fill="") +
scale_fill_manual(values=c('#f5200c','#fc9288'), labels = c("anterior", "posterior"))+
scale_x_discrete(limits=c("inside", "outside"), labels=c("inside window", "outside window"))+
scale_y_continuous(expand = c(0, 0), limits = c(0,0.3), breaks=seq(0,0.3,.05))+
geom_hline(yintercept=0)+
ggtitle("parahippocampal
supercluster")+
theme(plot.title = element_text(size=22, face="bold", hjust=0.5), legend.position="bottom")
#dev.off()
#ggsave("parahippocampal.png", path="~/Downloads", width = 4.5, height = 4.5)